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We investigate both continuous (second-order) and discontinuous (first-order) transitions to 
macroscopic synchronization within a single class of discrete, stochastic (globally) phase-coupled 
oscillators. We provide analytical and numerical evidence that the continuity of the transition de- 
pends on the coupling coefficients and, in some nonuniform populations, on the degree of quenched 
disorder. Hence, in a relatively simple setting this class of models exhibits the qualitative behaviors 
characteristic of a variety of considerably more complicated models. In addition, we study the mi- 
croscopic basis of synchronization above threshold and detail the counterintuitive subtleties relating 
measurements of time averaged frequencies and mean field oscillations. Most notably, we observe 
a state of suprathreshold partial synchronization in which time-averaged frequency measurements 
from individual oscillators do not correspond to the frequency of macroscopic oscillations observed 
in the population. 

PACS numbers: 64.60.Ht, 05.45.Xt, 89.75.-k 
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I. INTRODUCTION 

A great number of physical systems consist of indi- 
vidual entities with periodic, or nearly periodic, dynam- 
ics. Ranging from collections of chemical consitutents to 
groups of social entities - for example, applauding indi- 
viduals whose clapping is repetitive - these systems serve 
as a battleground of sorts for the competition between the 
dynamics of individual constituents and the large scale 
cooperation favored in many cases by the nature of their 
mutual interactions. Owing to the ubiquity and certainly, 
in part, to the dramatic nature of the emergent synchro- 
nized behavior in such naturally oscillating settings, the 
subject has been intensely studied in the physics liter- 
ature for several decades, with the Kuramoto oscillator 
and its kin serving asprototypical models on which many 
studies are based [3, [2, S 0| • 

Recently, a simple class of models of macroscopic syn- 
chronization have provided a number of additional in- 
sights into the large-scale phenomena ocurring in noisy 
discrete coupled oscillators, including detailed characteri- 
zations of both the universal critical behavior of the con- 
tinous phase transition 0, Q as well as the effects of 
spatial disorder in such populations Q. While retaining 
many qualitative characteristics of more complex mod- 
els, the discrete oscillators remain sufficiently simple to 
provide results unattainable in most of the paradigmatic 
settings. 

In this paper, we generalize our class of stochastic, 
discrete oscillator models and detail its use in a variety 
of new contexts. By generalizing the form of the inter- 
oscillator coupling, we show that our class of mean field 
models encompasses oscillators which can undergo either 
supercritical or subcritical Hopf bifurcations, depending 



on the microscopic specifics of the coupling. In addition, 
we study dichotomously disordered populations of oscil- 
lators and show that the bifurcation can be either super- 
critical or subcritical depending on the degree of disor- 
der in the population. Such behaviors are reminiscent of 
a number of si gnifi cantly more complex oscillator mod- 
els including Daido's generalized 
Kuramoto oscillators [151 ] , where either disoder or micro- 
scopic coupling specifics can alter the nature (continuous 
or discontinuous) of the transition. However, our model 
provides a far simpler setting for observing both conti- 
nous and discontinous transitions to synchronization. 

In addition, we study the microscopic underpinnings 
of synchronization above threshold. In particular, we 
look at time-averaged frequency and its relationship to 
phase synchronization above threshold (which turns out 
to be interestingly counterintuitive). Again, we do this 
for a specific model from our general class (for both sin- 
gle and dichotomously disordered populations), but we 
expect the results to hold for our entire class of models 
undergoing a Hopf bifurcation. This is somewhat similar 
to the partial synchronization seen in other models [l8|, 
but again, our model is simpler and, perhaps, more trans- 
parent. 

In Sec. [H] we present our model and highlight its es- 
sential parameters. Here we note that for systems of 
oscillators with identical transition rates between states 
the control parameter for the phase transition is the cou- 
pling strength among oscillators; when the array includes 
oscillators with different transition rates, the degree of 
disorder is also a control parameter. In Sec. IIIII we show 
that depending on the values of microscopic parameters, 
this model can exhibit both subcritical (or first-order) 
and supercritical (or second-order) phase transitions, as 
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a function of the coupling strength and also as a function 
of the degree of disorder. Section|IV]dcals with the micro- 
scopic underpinnings of the synchronization phenomenon 
and the connection between phase synchronization and 
frequency cntrainmcnt in our system. Section fVl presents 
a summary and further discussion of our results. 



II. THE MODEL 

In this section we present our model in some detail, re- 
peating some of our early presentations H H, 0] because 
of important (albeit simple) generalizations of the model. 
We begin by considering a stochastic three-state model 
governed by transition rates g (see Fig. [T|), where each 
state may be interpreted as a discrete phase @, H, 0|- 
Because the transitions among states are unidirectional 
and do not conform to deterministic rate laws, the model 
retains a qualitative link with a noisy phase oscillator. 
The linear evolution equation of a single oscillator is 
dP(t)/dt = MP[t), where the components Pj(t) of the 
column vector P(t) = (Pi(t) P 2 (t) P 3 (t)) T (T denotes 
the transpose) are the probabilities of being in states i 
at time t, where i — 1, 2, 3 and 
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(1) 



The system reaches a steady state for P* = P 2 * = P 3 * = 
1/3. The oscillator's periodicity, as contained in the 
timescale of the cycle i = 1 — » 2 — > 3 — » 1..., is de- 
termined by g; that is, the time evolution of our sim- 
ple model qualitatively resembles that of the discretized 
phase of a generic noisy oscillator with the intrinsic eigen- 
frequency set by the value of g. 

To study interacting arrays of these oscillators, we cou- 
ple individual units by allowing the transition rates of 
each unit to depend on the states of the units to which it 
is connected. Specifically, for N identical units we choose 
the transition rate of a unit v from state i to state i + 1 



as 



gi = g exp 



a(UN l+1 + VN^ + WN t ) 



(2) 



where i = 1,2,3 and i + 1 = 1 when i = 3, a is the 
coupling parameter, g is the transition rate parameter, n 




FIG. 1: Fhree-state unit with transition rates g. 



is the number of oscillators to which unit v is coupled, 
and iVfc is the number of units among the n that are in 
state k. We introduce the real constants U, V, and W 
to encompass in a general way our two previous coupling 
functions [1, d, 0j- Each unit may thus transition to 
the state ahead or remain in its current state, and the 
propensity for such a change depends on the states of 
its nearest neighbors. In our earlier works we considered 
the globally coupled system, n = N — 1, and also nearest 
neighbor coupling in square, cubic, or hypercubic arrays, 
n = 2d (d = dimensionality). Here we focus on globally 
coupled arrays. 

For a population of N — > oo identical units in the mean 
field (globally coupled) version of this model we can re- 
place Nk/N with the probability Pfc, thereby arriving 
at a nonlinear equation for the mean field probability, 
dP{t)/dt = M[P(t)]P(t), with 



-5i g 3 
M[P{t)\ = | 9l -g 2 

52 ~93, 



(3) 



Normalization allows us to eliminate Pj(t) and obtain a 
closed set of equations for Pi(t) and P2(t)- We can then 
linearize about the fixed point (P*,P 2 *) = (1/3,1/3), 
yielding a Jacobian A(a,g,U,V,W) with a set of com- 
plex conjugate eigenvalues which determine the stability 
of this asynchronous state. Specifically, we find that 



A± = C{ - 9 + 3aA uw 

± iV3(3 + a(U + W - 2V))) , 



(4) 



where C = ge a ( u+v+w ^ 3 /6 is a nonzero constant for all 
finite U, V, and W and we introduce the abbreviation 
A m „ = m — n. The eigenvalues cross the imaginary axis 
at a c = 3/Auwi yielding 



with 



X* ± = ±iu(U,V,W) 



w{U, V, W) = g ^(U+V+W)/A uw 



Auv 
A 



uw 



(5) 



(6) 



For A uw ^ and w(U,V,W) ^ (that is, A uv ^ 0), 
a c represents a Hopf bifurcation point, indicating the 
emergence of macroscopic oscillations indicative of syn- 
chronization. Furthermore, we require that Auw > 
to ensure the bifurcation happens at a positive value 
of a. We note that in previous studies we have used 
(U, V, W) = (1, -1, 0), yielding a c = 3andw = 2g^/2, 0], 
and (U,V,W) = (1,0,-1), yielding a c = 1.5 and u> = 
gV3/2 (HQ. In addition, we stress that while a range 
of models may prove useful for exploring the phase tran- 
sition behavior near threshold (see, for example, @, @]), 
only models with W — provide physically appealing 
characteristics far above threshold (see, for example, 0]). 
Specifically, only for W — does the frequency of a per- 
fectly synchronized set of oscillators maintain a nonzero 
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FIG. 2: (Color online) In a single population of globally cou- 
pled oscillators, two physically distinct Hopf transitions can 
be observed depending on the choices of U, V, and W. The 
top panel represents (U, V, W) = (1, — 2, 0) and clearly shows 
characteristics of a continous transition, including hysteresis. 
Squares represent solutions starting from ordered (mostly syn- 
chronized) initial conditions, while circles represent solutions 
starting from disordered (random) initial conditions. The 
bottom panel represents (U,V,W) — (2,-1,0) and displays 
a continous transition with critical exponent /3 given by the 
classical value 1/2. The inset shows a log-log plot near the 
critical point. For comparison, a dashed line with a slope of 
1/2 is shown along with the order parameter curve (solid line) 
to verify this scaling relation. 



finite value (g). Below we explore in more detail the na- 
ture of the Hopf bifurcation associated with the class of 
models described by the permitted values (U, V, W). 

In addition to the single population case, we con- 
sider globally coupled arrays of oscillators that can 
have one of Af < N different transition rate param- 
eters, g = 7„, u — l,...,Af. As detailed in Q, 
the probability vector is now 3A/"-dimensional, P{t) = 



(-Pl,7l -^2,71 f*3,7i 



Pi ,7Af ^2,7^ ^3,7aa) T j and the added 



subscript on the components of P(t) keeps track of the 
transition rate parameter. Explicitly, the component 
is the probability that a unit with transition rate 
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parameter g — 7 U is in state i. The evolution of the 
probability vector is given by the set of coupled nonlinear 
differential equations dP(t)/dt = Mj^[P(t)]P(t), with 
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M u [P{t)] = 



7i 
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72 



V o 



... 
... 

Mr* 



(7) 



Here 




33 (7u) 

92 (7u) ~53 {lu). 



(8) 



;(7u) = 7«exp 



M 
fe=i 



*(UP l+1 . lk +VP t - lnk +WP ink ) 



(9) 



The function <f("fk) is the fraction of units which have a 
transition rate parameter g = lk- 

Because it closely appeals to physical intuition [7| for 
oscillators far above threshold, we limit ourselves to the 
case (U,V,W) = (1,-1,0) for dichotomously disordered 
oscillators. We further limit our focus here to uniform 
distributions f{fyk) = 1/A/", but note that relaxing this 
constraint has been shown to preserve the qualitative fea- 
tures of the model Q. For uniform distributions and 
(U,V,W) = (1,-1,0), probability normalization again 
allows us to reduce this to a system of 2Af coupled ordi- 
nary differential equations. We can then linearize about 
the disordered state P(t) = (1/3 1/3 ... 1/3) T and ar- 
rive at a 2AT x 2Af Jacobian parameterized by a collec- 
tion of TV transition rate parameters 7$ and a coupling 
strength a. 

While it has been shown that the qualitative essence 
of the model remains similar for J\f = 2, 3, 4 and even for 
completely disordered populations Q , we focus here only 
on the simple dichotomously disordered case, N — 2. As 
shown in the four eigenvalues (A+, \* + , A_, AL) of the 
corresponding Jacobian are given by 



ReA± 

7i +72 
ImA± 

7i +72 
where 

and 



— [a — 6 ± B(a, /i) cos (C(a, /i))] 
8 

' f V3(a + 2) ±B(a,n) sin (C(a,/x)) 



(10) 



V2 [a 4 - 6aV + 3/(a 2 + 3)] 
'-\/3(a 2 - (a + 3)/i 2 )~ 



1/4 



1 

— tan 
2 



(11) 



a 2 + 3(a - l)n 2 



2|7i -72I 



(71+72) 



(12) 



Aside from an overall factor (71 +72), Eqs. (fT0|) depend 
only on the relative width variable /1, and therefore the 
critical coupling a c , that is, the value of a at which Re 
A + = 0, depends only on /1. As Re A_ does not vanish 
for any a, a c corresponds to a Hopf bifurcation, and our 
J\f = 2 model exhibits macroscopic oscillations indica- 
tive of large-scale cooperation. We note that a c increases 
with increasing /1, indicating that a stronger coupling is 
necessary to overcome increasingly different values of ji 
and 72. 
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In what follows, we make use of the synchrony order 
paramter r to characterize the emergence of phase syn- 
chrony [TiT ]. This parameter is defined as 



1 N 

r=(R), i^|£< 

3=1 



I 



(13) 



Here <fi is a discrete phase, taken to be 2w(k — l)/3 for 
state A: € {1,2,3} at site j. The brackets represent an 
average over time in the steady state and an average over 
all independent trials. Therefore, r serves as a measure 
of phase synchronization. 



III. CONTINUOUS AND DISCONTINUOUS 
TRANSITIONS TO SYNCHRONY 

In the mean field limit, the order of the phase transi- 
tion to synchrony is closely tied to the nature of the Hopf 
bifurcation. Specifically, a subcritical Hopf bifurcation 
corresponds to a discontinuous (sometimes called first- 
order) phase transition, while a supercritical Hopf bifur- 
cation indicates a continuous (second-order) transition. 
As such, we place special emphasis in this manuscript on 
the sign of l±, the first Lyapunov coefficient, which pro- 
vides information on the nature of the Hopf bifurcation 
and, by extension, on the order of the phase transition. 

In general, l\ can be calculated using the projection 
technique given in 16], which relies on a multivariate 
Taylor expansion of the vector field describing the dy- 
namics in question about an equilibrium point. For a 
general rt-dimcnsional dynamical system x — f(x,e) with 
an equilibrium point x — x H undergoing a Hopf bifurca- 
tion at parameter value e — e H , l\ is given by (l6j 
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h = ^-Re«P, C(q, q, q)) - 2<p, B(q, A^B^ q))) 



(14) 



+ {p,B(q : (2luI-A)- 1 B(q,q)))), 
where (., .) is the typical complex scalar product, / is the 
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FIG. 3: The first Lyapunov coefficient h is shown for Hopf 
bifurcations taking place at e H = (a c (n),^)- The bifurcation 
can be either subcritical or supercritical depending on the 
relative width variable. 



identity matrix, and p and q are right and left eigenvec- 
tors of the Jacobian A = ^- 1 „ given by 

ox \x—x H ° J 



Aq = itoq 
A T p = —icup. 



(15) 



Furthermore, p is chosen so that (p,q) — 1, and B(u,v) 
and C(u, v, w) are multilinear, n-dimensional vector func- 
tions corresponding to the lowest order nonlinear coeffi- 
cients in the Taylor expansion of the vector field. That 



B(u,v) 



E 

j,k=i 



UjVk, 



ip—x 1 



C{u,v,w)= 



(16) 



UjVkWt, 



tp—x 1 



with x H indicating the equilibrium point of the vector 
field around which we expand and e H the bifurcation 
parameter, e, evaluated at the bifurcation point. 



A. Continuous and Discontinuous Transitions in a 
Single Population of Identical Oscillators 

For the case of a single population of oscillators de- 
scribed by Eqs. ^ and ([3]), Zi can be analytically cal- 
culated using the technique outlined above. Specifically, 
we set g — 1 (without loss of generality) and consider the 




FIG. 4: (Color online) A subcritical Hopf bifurcation occurs 
for jj, — 3/4. Squares represent solutions starting from or- 
dered (mostly synchronized) initial conditions, while circles 
represent solutions starting from disordered (random) initial 
conditions. Dark (blue) points correspond to population one, 
71 = 2.5, and light (pink) points to population two, 72 = 5.5. 
The transition is clearly discontinuous as a crosses a c ~ 3.55. 
In addition, a region of multistability and corresponding hys- 
teresis exists just below threshold. 
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equilibrium point P = (1/3, 1/3) at e H = a c and find q 
and p to be 



(17) 




V3 2i 



+ V3' Zi + ^/Z)' 

independent of Z7, V, and Then, calculating the mul- 
tivariable functions B(u, v) and C(u,v,w) with Eq. (116|) 
and using w as defined in Eq. ([6]) along with Eqs. (fM)) 
and p7[) . we find after simplification that 



Zi = - 



9V3(C7 + V - 2W) 



4A 



u w 



(18) 



As a result, the nature of the Hopf bifurcation depends 
on the choices U, V, and W. Specifically, if we assume 
U > V, we have 



h < for IT < 
Zi > for IT > 



u + v 
u + v 



(19) 



A similar result holds for U < V, but we shall here re- 
strict ourselves to the intuitively reasonable models posit- 
ing U > and V < 0; that is, the oscillators one state 
ahead of the one in question can only increase (or not 
affect) the transition rate and those behind can only de- 
crease (or not affect) the transition rate. To verify these 
predictions, we show numerical solutions to the mean 
field equations in Fig. [51 the top panel represents an ex- 
ample in the subcritical regime ((Z7, V, W) = (1,-2,0)) 
while the bottom panel shows an example in the su- 
percritical regime ((U,V,W) — (2,-1,0)). A clear dis- 
tinction can be made in the neighborhood of the crit- 
ical point. We also note that the continuous transi- 
tion is characterized by the classical mean field exponent 
0=1/2. 

We further observe that the choice (U, V, IT) = 
(1,0,-1) leads to h = -27^3/4 « -11.69, indicat- 
ing a supercritical Hopf bifurcation and rendering the 
model applicable to studies of continuous phase transi- 
tions Q. With universality in mind, we stress that 
any choice of parameters (U, V, W) yielding a supercriti- 
cal bifurcation should show similar critical behavior. On 
the other hand, the choice (U,V,W) = (1,-1,0), while 
physically appealing above threshold, falls at a singular 
point separating the subcritical and supercritical cases 
(Zi = 0). The flexibility inherent in the choice of coeffi- 
cients X, Y, and Z speaks to the richness of our generic 
three-state oscillator and highlights its utility in studying 
synchronization in both supercritical (see [1,[6[) and sub- 
critical regimes. We proceed to study the model in the 
presence of dichotomous disorder and show that, for a 
given choice (£/, V, W), the level of disorder can also alter 
the nature of the Hopf bifurcation and hence the order of 



the phase transition. We select the physically appealing 
choice (1, —1, 0) and, while the behavior near the critical 
point may depend in some sense on this choice, we stress 
that our overarching goal remains unhindered. That is, 
we are able to provide an example indicating that dis- 
order alone can affect the nature of the transition. The 
ubiquity of this phenomenon across the entire range of 
models remains an open question for future work, though 
we note that similar results are observed for all parameter 
choices mentioned in this paper. 



B. Continuous and Discontinuous Transitions in a 
Dichotomously Disordered Population 

Interestingly, the dichotomously disordered system 
corresponding to Eqs. 0, ©, and © with N = 2 
and (U, V, W) = (1, —1, 0) can undergo either a subcrit- 
ical or supercritical bifurcation depending on the value 
of /i characterizing the individual transition rates. The 
transition to synchrony occurs at a single value of the 
coupling a c (n) dependent on the relative width param- 
eter Q- As such, a and fi are not truly independent 
parameters, and we can in principle eliminate a and con- 
sider fi to be the bifurcation parameter of interest. Then, 
using the machinery of Eqs. ([14)). (fl~5)) , and (fT6j). it is a 
straightforward but tedious exercise to numerically eval- 
uate the first Lyapunov coefficient h{n) corresponding 




FIG. 5: (Color online) A supercritical Hopf bifurcation occurs 
for n = 7/4. Squares represent solutions starting from an 
ordered (mostly synchronized) initial condition, while circles 
represent solutions starting from a disordered (random) initial 
condition (repetition with other ordered and disordered initial 
conditions leads to essentially the same results). Dark (blue) 
points correspond to population one, 71 = 0.25, and light 
(pink) points to population two, 72 = 3.75. The transition 
is clearly continuous as a crosses a c ta 5.44, and there is a 
noticeable absence of hysteresis. Dotted line is drawn to guide 
the eye. The inset indicates that, in the neighborhood of the 
critical point, the order parameter follows power law behavior 
with the correct mean field critical exponent (/3 = 1/2). 
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to the Hopf bifurcation occurring at (/x, a(/i)). As shown 
in Fig. [31 the sign of varies depending on the relative 
width parameter (which in turn determines the critical 
coupling a c ). Hence, the phase transition to synchrony 
can appear continuous or discontinuous depending on the 
relative difference between the transition rate parameters 
in the two populations. 

To verify these predictions, we solve the mean field 
equations numerically in both the subcritical (/i = 3/4) 
and supercritical (/i = 7/4) regimes. In the former case, 
we consider the case 71 = 2.5, 72 = 5.5. Figure [4] clearly 
indicates that the transition to synchrony is marked by 
a discontinuous change in the order parameter r as a 
eclipses a c ~ 3.55. In addition, a small region of marked 
hysteresis appears just below threshold. Remarkably, 
this indicates that a stable disordered solution coexists 
with a stable, synchronized solution (the stable limit cy- 
cle) just before threshold. 

By contrast, the case [i = 7/4 corresponds to a su- 
percritical Hopf bifurcation reminiscent of a continuous 
phase transition. As shown in Fig [5l the transition is 
characterized by a continously increasing order parame- 
ter; no hysteresis is evident. We note also that the order 
parameter displays a power law increase near the onset 
of the bifurcation marked by the mean field critical ex- 
ponent (3 — 1/2. This is expected both from the Hopf 
bifurcation theorem, which prescribes the (a — ac) 1 ^ 2 de- 
pendence of the limit cycle radius (closely related to r, 
the order parameter) near the onset of synchrony, and 
also because of the analogy with phase transitions in an 
infinite-dimensional space (see 0,11). 

Interestingly, these results indicate that the degree of 
spatial disorder may fundamentally alter the nature of 
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FIG. 6: (Color online) Each plot shows a histogram of time- 
averaged frequencies (in the steady state), where the vertical 
axis represents the number of units (out of TV = 3500 total 
units) having the frequency ui. The power spectrum of -Pi, 71 
overlays each histogram. The top panel is below synchroniza- 
tion threshold (a = 2.65), while the middle (a = 3.05) and 
lower panels (a = 3.45) are both above threshold. 



the phase transition to synchrony. In both the subcriti- 
cal and supercritical cases, synchronization is marked by 
the destabilization of the non-synchronous state at a sin- 
gle value of a c , giving rise to emergent oscillations in a 
macroscopic variable [for example, P\ ni {i)\. Hence, both 
cases retain the qualitative features of synchronization in 
disordered populations discussed in our previous work [7j ; 
however, the details of the onset of such cooperation dis- 
tinguish the two cases. 



IV. MICROSCOPIC UNDERPINNINGS OF 
SYNCHRONIZATION 

Having detailed two distinct mechanisms by which syn- 
chronization might arise, we now explore in detail the 
microscopic subtleties underlying synchronization above 
threshold. As detailed in [f| [y, an d mentioned above, 
synchronization occurs in the mean field limit via the 
destabilization of a nonsynchronous fixed point. Specifi- 
cally, a single pair of complex conjugate eigenvalues cor- 
responding to the linearized fixed point cross the imag- 
inary axis at <z c , giving rise to stable oscillations in the 
macroscopic variables characterizing the system [in our 
case, the components of P(t)]. While the onset of this 
behavior is dependent on the choice (U, V, W) and also 
the magnitude of disorder within the system (Sec. IIII|I . 
the qualitative features of the synchronized state re- 
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FIG. 7: (Color online) Each plot shows a histogram of time- 
averaged frequencies (in the steady state), where the vertical 
axis represents the number of units (out of TV = 3500 total 
units) having the frequency Q. Population one, characterized 
by 71 = 2.5, is represented by the dark (blue) histogram, while 
population two, characterized by 72 = 5.5, is represented by 
the light (pink) one. Power spectra of Pi, 71 (dark, blue) and 
Pi, 72 (light, pink) overlay the histograms. The top panel is 
below synchronization threshold (a — 3.20), while the middle 
(a = 3.60) and lower panels (a — 3.86) are just above thresh- 
old. For aesthetic purposes, the horizontal range is relatively 
shifted in the three panels. 
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main identical above threshold in both the subcritical 
and supercritical case. Hence, we limit our attention 
to several illustrative cases, but note that our results 
hold also for the supercritical case (and in fact then en- 
tire range of /x). Specifically, in what follows, we take 
(U,V,W) = (1,-1,0) and consider a single population 
as well as a dichotomously disordered population with 
fi = 3/4. 

In particular, the threshold a c is marked by the on- 
set of coherent temporal oscillations in the components 
of P(t). We characterize the microscopic underpinnings 
of these oscillations by considering Qi , the time-averaged 
frequency of oscillator i in the steady state. We per- 
form simulations on globally coupled lattices of N = 3500 
units of a single population with 7 = 1 and also of a di- 
chotomously disordered population with 71 = 2.5 and 
72 = 5.5. As shown in Figs. [5] and [7J the distribution of 
frequencies uJi clusters around the values prescribed by 7 
(or 71 and 72 for populations one and two, respectively) 
far below threshold (top panels). Specifically, for a de- 
terministic oscillator with transition rate 7, UJi is given 
by 27T7/3; when 7=1 (or 71 = 2.5 and 72 = 5.5), 
this gives the central peak of the histogram for the rele- 
vant population. We compare these histograms with the 

/- - \l/2 

power spectra ( Pi i7i (w)P* )7 . (u)) ) , where Pi )7i (w) is 

the Fourier transform of Pi 7i (t). As threshold is eclipsed 
(middle panels), a peak arises in the power spectrum of 
the macroscopic variables Pi, 7i , though the frequency of 
this peak does not correspond with the individual cJj's of 
oscillators constituting the population. In the dichoto- 
mous case, this peak only roughly corresponds with the 
time averaged frequencies from population two and com- 
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FIG. 8: (Color online) The central figure shows a histogram of 
waiting times Ti step for a sampling of N — 60 units once the 
steady state has been reached. The vertical black line indi- 
cates the waiting time corresponding to the peak in the power 
spectrum of Pi(i) (that is, the waiting time corresponding 
to the frequency of the macroscopic oscillations). The inset 
shows a similar histogram for a single unit, a — 3.05 (above 
threshold), 7 = 1 in all plots. 




FIG. 9: (Color online) The central figure shows a histogram 
of waiting times Ti 3tep for a sampling of N = 30 units from 
each population once the steady state has been reached. Dark 
(blue) represents population one, while light (red) represents 
population 2. The vertical black line indicates the waiting 
time corresponding to the peak in the power spectrum of Pi (t) 
(that is, the waiting time corresponding to the frequency of 
the macroscopic oscillations). The insets show similar his- 
tograms for single units; the top histogram is for a unit from 
population one, the bottom from population two. a — 3.60 
(above threshold), 71 = 2.5, and 72 = 5.5 in all plots. 



pletely exceeds even the maximum UJi characterizing pop- 
ulation one. As a is further increased, the descrepency 
between the time-averaged frequency histograms and the 
macroscopic oscillation frequencies decreases. In addi- 
tion, in the disordered case, the histograms for the two 
populations become increasingly narrow and closer to one 
another (bottom panel). We note that as a becomes 
tremendously large, the histograms become extremely 
narrow and begin to overlap at a frequency determined by 
the frequency of the macroscopic oscillations, as expected 
(indicative of perfect synchronization). Nonetheless, the 
behavior for finite, intermediate a is rather counterintu- 
itive and points to a rich microscopic dynamics underly- 
ing the cooperative behavior. 

To further explore these trends, we consider the 
stochastic variable T± step , the waiting time in a sin- 
gle state for an individual oscillator. T\ s t ep represents 
the time the oscillator spends in a single state i be- 
fore transitioning to the subsequent state i + 1. For 
computational efficiency, we record T\ step for a repre- 
sentative subpopulation of 60 units (30 units from each 
population, one and two, in the disordered case). Fig- 
ures [8] and [9] show histograms of the variable T\ step taken 
over this representative subpopulation once steady state 
was reached. Clearly, all relevant subpopulations con- 
sist of oscillators whose steps most often correspond to 
the frequency of the macroscopic oscillation (shown by 
the solid vertical line). That is, the peak of the his- 
tograms occur at a value T comensurate with the fre- 
quency peak in the power spectrum of the components 
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of P(t). However, Fig. [8] shows that the distribution of 
T\ step is bimodal, with a significant peak occuring at 
T\ step ~ 2.2 which downward biases the time-averaged 
frequencies oji of individual units. We note that as cou- 
pling a increases significantly above threshold, the distri- 
bution becomes unimodal with a peak at T\ step corre- 
sponding to the frequency of macroscopic oscillation. In 
the disordered case, only population one, characterized 
by significantly lower time-averaged tD's, shows a bimodal 
distribution with a significant peak at T\ s tep ~ 0.45. In 
fact, these long waiting times, while not the dominant 
macroscopic behavior, pervade the microscopic dynam- 
ics in such a way that the time-averaged frequencies be- 
come downward biased and no longer accurately repre- 
sent the macroscopic dynamics. Interestingly, population 
two has become sufficiently synchronized that the second 
peak is effectively nonexistent, and thus the frequencies 
overlap more closely with the macroscopic "mean field" 
frequency. The right insets of Figs. [5] and [5] show his- 
tograms for single units chosen from the populations (or 
subpopulations). Again, the unit chosen from the single 
population case shows a bimodal distribution with sig- 
nificant "anamolous" peaks near T\ step ~ 2.2. In the 
disordered case, the unit from population one shows a 
bimodal waiting time distribution characterized by oc- 
casional waiting times in the neighborhood of T rs 0.45 
in addition to those corresponding to the macroscopic 
oscillations. Finally, in Fig. [10] we show the time evo- 
lution of the subpopulations along with the macroscopic 
variable Pi )7j for each population. At any given time, 
the majority of oscillators in each population is synchro- 
nized, leading to the smooth oscillations of the macro- 
scopic variable. However, isolated single units are prone 
to long waiting times, particularly in the less synchro- 
nized population (population one, left panel, in this ex- 
ample). These anamolously long waiting times, which 
serve to bias the time averaged frequencies Ui of each in- 
dividual unit, nevertheless do not substantially disrupt 
the macroscopic oscillations, largely because the occu- 
rance of coincident long waits is fairly uncommon. That 
is, the long waiting times do not appear in any signif- 
icantly correlated way among individual constituents of 
the population. 



V. DISCUSSION 

We have shown that a class of simple, discrete models 
of stochastic phase coupled oscillators can undergo either 
a subcritical or supercritical bifurcation to macroscopic 
synchrony, depending on the chosen form of the micro- 
scopic coupling. As such, the different instances of the 
model can be used to study either continuous phase tran- 
sitions 0, Q or discontinuous transitions exhibiting hys- 
teresis, a characterstic seen in detailed theoretical mod- 
els of, e.g., coupled Josephson junctions 0] but only ob- 
served in significantly more complex coupled oscillator 
models We stress that uni- 






FIG. 10: (Color online) The top panels show the evolution 
of a representative sub-system (N — 30 units of each popu- 
lation). Dark, medium, and light represent states 1,2, and 3, 
respectively. The bottom panels show the macroscopic vari- 
able Pi, 7i for each population. The left panels correspond to 
population one (71 = 2.5), while the right panels correspond 
to population two (72 = 5.5). a — 3.60 (above threshold) for 
all plots. 



versality suggests that all models in this class exhibiting 
continous phase transitions should show similar behavior 
near the criticalpoint, and this served as the basis of our 
eariler studies [a, @]. Nevertheless, it is remarkable that 
minor modifications in microscopic coupling can alter the 
nature of the bifurcation in such a fundamental way. 

In addition, we have shown that in dichotomously dis- 
ordered populations, both subcritical and supercritical 
Hopf bifurcations can occur, and the distinction is com- 
pletely determined by the relative width /x characterizing 
the transition rate disorder between the two populations. 
While the qualitative features of the transitions within 
each class (subcritical and supercritical) appear identical, 
the distinction between classes points to fundamentally 
different mechanisms underlying the initial emergence of 
phase synchronization. In particular, it is striking that 
the level of disorder within a population, as measured by 
/i, can significantly alter the behavior near the critical 
point (though we stress that behavior even moderately 
above threshold is qualitatively indistinguishable). 

Finally, we have studied the microscopic basis of phase 
synchronization above threshold. It is initially counter- 
intuitive that phase synchronization, defined in terms 
of the Hopf bifurcation and temporal oscillations in the 
macroscopic variable P(t) (and measured in the order 
parameter r), is not contingent upon the existence of 
overlapping distributions of oJi. That is, our results re- 
garding the discrete oscillator model highlight the com- 
plexity of microscopic dynamics underlying macroscopic 
cooperation and point to a potentially misleading sub- 
tlety. Whereas phase synchronization is often considered 
a stronger condition than frequency entrainment - de- 
fined using an order parameter built upon the notion 
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that a fraction of units display identical time-averaged 
frequencies in the oscillator population - we here re- 
port subtle microscopic features which distinguish the 
two without establishing a clear hierarchy. For exam- 
ple, Hong et al. [Tjj show that for disordered popula- 
tions of Kuramoto oscillators, the lower critical dimen- 
sion for frequency entrainment is lower than that for 
phase synchronization in locally coupled oscillators, in- 
dicating the relative ease with which frequency entrain- 
ment is achieved. They note that the two transitions 
coincide in the case of globally coupled units. Contrast 
that with our dichotomously disordered population, for 
which phase synchronization occurs without any overlap 
in the frequency distributions: that is, no oscillator from 
population one has the same frequency as any oscilla- 
tor from population two. While a direct comparison is 
not plausible owing to the specific differences between 
models and order parameters, we stress that any order 
parameter related to time-averaged measurements of fre- 
quencies would, for our model, be misleading and provide 
potentially counterintuitive results. The emergence of a 
nonzero r, which measures phase synchronization, cor- 
responds with the loss of stability of the asynchronous 
fixed point (the Hopf bifurcation). This does not guar- 
antee similar distributions of time-averaged frequencies 
in the two populations; in fact, we can readily see that 
synchronization occurs while the frequency distributions 
are entirely distinct. Furthermore, the frequency of the 
macroscopic oscillations of the mean field does not always 
coincide with the time-averaged frequencies of the oscilla- 
tors constituting the population (or any subpopulation). 
Only when coupling is sufficiently large to substantially 
reduce the anamolously long waiting times which bias uJi 
will the frequency distributions begin to overlap one an- 
other and coincide with the frequency of the mean field 
oscillations. Because these long waiting times appear 
more readily in the population with the smaller 7$, the 
time-averaged frequencies of the two populations are dis- 
proportionately affected, meaning that the populations 
will appear to behave quite differently in terms of aver- 
age frequency. This in fact underlies the stark differences 
in the degree of synchronization between two populations 
as measured by the order parameter r, and provides an 
intuitive description capable of explaining this discrep- 



ancy. Our previous results show that completely disor- 
dered populations show qualitative similarities with the 
dichotomously disordered case Q; hence, we are led to 
cautiously speculate that wholly disordered populations 
are also characterized by waiting times T\ step distributed 
with long tails, and hence time-averaged frequencies be- 
come downwardly biased, meaning that the order param- 
eter for frequency entrainment, in the typical sense, will 
not accurately reflect the macroscopic cooperation. Fur- 
ther studies along these lines are currently in progress. 

Finally, the results of this work raise the following 
question: how dependent is the above phenomenon on 
the choice of a discrete phase model? Would simi- 
larly counterintuitve results arise in continous phase os- 
cillators? In fact, a recent study by Rosenblum and 
Pikovsky [r| suggests that a similar (though not iden- 
tical) state of partial synchronization arises in continous 
oscillators coupled in a highly nonlinear fashion. Specifi- 
cally, they find that in globally coupled oscillators, phases 
exist in which certain subpopulations are characterized 
by time-averaged frequencies which are not commen- 
surate with the oscillations of the mean field, that is, 
they are not locked with the macroscopic oscillations in- 
duced in the population. While once again the differ- 
ences between the models make direct comparison diffi- 
cult, it is nevertheless clear that measurements of time- 
averaged frequencies provide potentially counterintuitive 
results, even in globally coupled arrays. In the case of 
our stochastic discrete oscillators, the behvavior is quite 
transparent once viewed in terms of T\ step , though it 
is not clear whether a similar mechanism underlies the 
phenomenon in the continous phase model. Uncovering 
the relationship between the superthrcshold phase in our 
model and that in the continous oscillator model of [HI 
remains an open question, but even the superficial sim- 
ilarities between the results motivate continued efforts 
along these lines. 
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